'''
Rough estimate of slinky sensitivity using Kingman, Rowland, and
Popescu (Am. J. Phys, 70(6), June 2002, pp 595-598) and estimates for
the magnetic properties of a K&J Magnetics RC44, courtesy of Michael Paul.

usage: python sensitivity.py [options]
options:
    -v                  be more verbose
    -h                  print this and exit
    -Z min,delta,max    range of Z to examine (cm)
'''

import sys
import traceback
import numpy
import matplotlib.pyplot as pyplot


mext = 1.65 # RC44 A/m**2 from Michael Paul


class SingleTurn(object):
    '''Voltage response of a single turn, of radius rm (meters),
    located a distance d below the magnet.  The calculation
    becomes inaccurate as rm approaches the radius of the magnet.'''

    def __init__(self, rm, m = mext):
        self.rm = rm
        self.m = m
        self.u0 = 1.256637061e-6 # N/A**2
        self.a = 3.0 * self.u0 * self.m / (2.0 * self.rm**2)

    def V(self, d, v = 1.0):
        zeta = d / self.rm
        return self.a * zeta * v / (1.0 + zeta**2)**2.5



class FullCoil(object):
    '''Voltage response to a magnet moving at unit velocity (1.0 m/s).
    Sensor has N turns and extends from z = 0.0 to z = -L.'''

    def __init__(self, rm = 0.0254, N = 4000, L = 0.0254, m = mext):
        self.rm = rm
        self.N = N
        self.L = L
        self.m = m
        self.turn = SingleTurn(self.rm, self.m)


    def V(self, d):
        '''Just sum the contributions of the individual turns.'''
        v = 0
        dz = self.L / self.N
        for nu in range(self.N):
            znu = float(nu) * self.L / self.N
            v += self.turn.V(d + znu, 1.0)
        return v
        


def main(argv=None):
    import getopt

    if argv is None:
        argv = sys.argv
    options = "hvZ:"
    verbose = 0
    Z = 1.0e-2 * numpy.arange(2.0, 0.0, -0.05)

    try:
        try:
            opts, args = getopt.getopt(argv[1:], options, ["help"])
            for opt, val in opts:
                if opt == '-h':
                    print __doc__
                    return 0
                elif opt == "-v":
                    verbose += 1
                elif opt == "-Z":
                    zs, zd, ze = [1.0e-2 * float(s) for s in val.split(',')]
                    Z = numpy.arange(zs, ze, zd)
                else:
                    raise RuntimeError("cannot handle option '%s'" % opt)

        except getopt.error, msg:
            print >> sys.stderr, __doc__

        v = []
        zl = []
        ring = FullCoil(rm = 0.0254, N = 4000, L = 0.0254, m = mext)
        for z in Z:
            v.append(ring.V(z))
            zl.append(100.0 * z) # convert to centimeters for plot
            if verbose > 1:
                print >> sys.stderr, "%11.4e  %10.3e" % (zl[-1], v[-1])

        pyplot.title("induced voltage for unit speed")
        pyplot.xlabel("cm")
        pyplot.ylabel("V")
        pyplot.plot(zl, v)
        pyplot.show()

    except RuntimeError, err:
        print >>sys.stderr, err
        return 2

    except:
        traceback.print_exc()
        return 1

if __name__ == "__main__":
    sys.exit(main())
    
    
